    ttime=ttime+dt2
	pNew=0.0d0
	Tnew=0.0d0
	Unew=0.0d0
	Enew=0.0d0
	rhonew=0.0d0
! viscosity 
	gradU=0.0d0
! ****************************************************************************************************************************** !
!	do k=1,nInternalFaces
!		DynVisc = DynamicViscosity(Ts(k), Rhos(k))
!		BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!		BulkViscousStress = BulkViscTotal * & ! BulkVisc * (V,Sf)
!				sum((U(1:3,meshFaceOwner(k))-U(1:3,meshFaceNeighbour(k)))*meshSf(1:3,k))
!		do i=1,3
!			delta=U(i,meshFaceOwner(k))-U(i,meshFaceNeighbour(k))
!			do j=1,3
!				vtmp=DynVisc * meshSf(j,k)*delta
!				do l=1,4
!					gradU(i,j,meshFaces(l,k))=gradU(i,j,meshFaces(l,k))+vtmp ! m2 * m/s
!					if(i.eq.j)then
!						gradU(i,j,meshFaces(l,k)) = gradU(i,j,meshFaces(l,k)) + BulkViscousStress
!					end if
!				end do
!			end do
!		end do
!	end do
!
!	do j=1,nBoundaryFields
!
!		if(btype(j).eq.3)then
!
!			do k=boundaryMeshStart(j)+1,boundaryMeshStart(j)+boundaryMeshSize(j)
!				DynVisc = DynamicViscosity(Ts(k), Rhos(k))
!				BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!				BulkViscousStress = 0d0
!
!				do i = 1,3
!					m = opposingOwnerFaceLabel(k)
!					delta = (4d0*U(i,meshFaceOwner(k)) - 0.5d0*(us(i,m)+usnew(i,m))) / 2d0
!					BulkViscousStress = BulkViscousStress + BulkViscTotal * delta * meshSf(i,k)
!				end do
!
!				do i=1,3
!					m = opposingOwnerFaceLabel(k)
!					delta = (4d0*U(i,meshFaceOwner(k)) - 0.5d0*(us(i,m)+usnew(i,m))) / 2d0
!					do m=1,3
!						vtmp = DynVisc * meshSf(m,k)*delta
!						do l=1,4
!!							iPnt = meshFaces(l,k)
!							gradU(i,m,meshFaces(l,k))=gradU(i,m,meshFaces(l,k))+vtmp ! m2 * m/s
!							if(i.eq.m)then
!								gradU(i,m,meshFaces(l,k)) = gradU(i,m,meshFaces(l,k)) + BulkViscousStress
!							end if
!						end do
!					end do
!				end do
!
!			end do
!		end if
!
!		if(btype(j).eq.4)then
!cycle
!			do k=boundaryMeshStart(j)+1,boundaryMeshStart(j)+boundaryMeshSize(j)
!				DynVisc = DynamicViscosity(Ts(k), Rhos(k))
!				BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!				BulkViscousStress = 0d0
!
!				do i = 1,3
!					m = opposingOwnerFaceLabel(k)
!					delta = ( 4d0*U(i,meshFaceOwner(k)) - 0.5d0*(us(i,m)+usnew(i,m) + 3d0*(us(i,k)+usnew(i,k)) ) ) / 2d0
!					BulkViscousStress = BulkViscousStress + BulkViscTotal * delta * meshSf(i,k)
!				end do
!
!				do i=1,3
!					m = opposingOwnerFaceLabel(k)
!					delta = (4d0*U(i,meshFaceOwner(k)) - 0.5d0*(us(i,m)+usnew(i,m) + 3d0*(us(i,k)+usnew(i,k)) ) ) / 2d0
!					do m=1,3
!						vtmp = DynVisc * meshSf(m,k)*delta
!						do l=1,4
!!							iPnt = meshFaces(l,k)
!							gradU(i,m,meshFaces(l,k))=gradU(i,m,meshFaces(l,k))+vtmp ! m2 * m/s
!							if(i.eq.m)then
!								gradU(i,m,meshFaces(l,k)) = gradU(i,m,meshFaces(l,k)) + BulkViscousStress
!							end if
!						end do
!					end do
!				end do
!
!			end do
!			cycle
!		end if
!
!!!		if (btype(j).ne.3.and.btype(j).ne.4.and.btype(j).ne.5) then ! INLET, OUTLET, PROC
!!		if (btype(j).ne.5) then ! INLET, OUTLET, PROC ! btype(j).ne.4.and.
!!!			if(btype(j).eq.1.or.btype(j).eq.2)then ! INLET, OUTLET
!!			if(btype(j).eq.1.or.btype(j).eq.2.or.btype(j).eq.3.or.btype(j).eq.4)then ! INLET, OUTLET
!!				factor=0.0d0 ! NOTE: factor=1.0d0 give mistake: wrong velocity on the boundary, VSemiletov
!!			else
!!				factor=1.0d0	! PROC boundary
!!			end if
!
!		if(btype(j).eq.1.or.btype(j).eq.2)then ! INLET, OUTLET
!			factor=0.0d0 ! NOTE: factor=1.0d0 give mistake: wrong velocity on the boundary, VSemiletov
!
!			do k=boundaryMeshStart(j)+1,boundaryMeshStart(j)+boundaryMeshSize(j)
!				DynVisc = DynamicViscosity(Ts(k), Rhos(k))
!				BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!				BulkViscousStress = 0d0
!				do i = 1,3
!					delta=factor*U(i,meshFaceOwner(k)) + (1.0d0-factor)*(U(i,meshFaceOwner(k)) - 0.5d0*(us(i,k)+usnew(i,k)))
!!					delta=factor*U(i,meshFaceOwner(k)) + (1.0d0-factor)*(U(i,meshFaceOwner(k)) - us(i,k))
!					BulkViscousStress = BulkViscousStress + BulkViscTotal * delta * meshSf(i,k) ! BulkVisc * (V,Sf)
!				end do
!
!
!				do i=1,3
!					delta=factor*U(i,meshFaceOwner(k)) + (1.0d0-factor)*(U(i,meshFaceOwner(k)) - 0.5d0*(us(i,k)+usnew(i,k))) ! Corrected BC, SK
!!					delta=factor*U(i,meshFaceOwner(k)) + (1.0d0-factor)*(U(i,meshFaceOwner(k)) - us(i,k)) ! Corrected BC, SK
!!					delta=factor*U(i,meshFaceOwner(k))
!					do m=1,3
!						vtmp = DynVisc * meshSf(m,k)*delta
!						do l=1,4
!							gradU(i,m,meshFaces(l,k))=gradU(i,m,meshFaces(l,k))+vtmp ! m2 * m/s
!							if(i.eq.m)then
!								gradU(i,m,meshFaces(l,k)) = gradU(i,m,meshFaces(l,k)) + BulkViscousStress
!							end if
!						end do
!					end do
!				end do
!			end do
!		end if
!
!!		if (btype(j).eq.4) then
!!			do k=boundaryMeshStart(j)+1,boundaryMeshStart(j)+boundaryMeshSize(j)
!!				i1=meshFaceOwner(k)	
!!				!	vtmp=dsqrt(U(3,i1)**2+U(3,i1)**2+U(3,i1)**2)
!!				!	localRe=Db(k)/2d0*vtmp/nu
!!
!!
!!				!	factor=rwall
!!				!*defFactor(localRe,vtmp,nu)
!!				do i=1,3
!!					! SK MODIFICATION FOR THE SECOND-ORDER VISCOUS WALL TREATMENT  
!!
!!						delta=factor*U(i,meshFaceOwner(k)) + 2.0 * (1.0d0-factor)*(U(i,meshFaceOwner(k)) - us(i,k))  ! Corrected BC, SK
!!!						delta=6.0d0*U(i,i1)-2.0d0*us(i,opposingOwnerFaceLabel(k))
!!
!!					!  ANOTHER POSSIBLE WAY(1):   =4.0d0*U(i,i1)-1.0d0*us(i,opposingOwnerFaceLabel(i))
!!
!!					! ORIGINAL MA'S FIRST ORDER: =factor*U(i,i1)
!!
!!					! END SK MODIFICATION 
!!					do m=1,3
!!						vtmp=meshSf(m,k)*delta
!!						do l=1,4
!!							gradU(i,m,meshFaces(l,k))=gradU(i,m,meshFaces(l,k))+vtmp ! m2 * m/s
!!						end do
!!					end do
!!				end do
!!			end do
!!		end if
!
!	end do
!!   periodic
!	do j=1,nBoundaryFields
!        if (btype(j).eq.5) then
!			ndeltap=boundaryMeshSize(j)/2
! 			do k=1+boundaryMeshStart(j),boundaryMeshStart(j)+ndeltap
!				DynVisc = DynamicViscosity(Ts(k), Rhos(k))
!				BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!				BulkViscousStress = BulkViscTotal * & ! BulkVisc * (V,Sf)
!						sum((U(1:3,meshFaceOwner(k))-U(1:3,meshFaceOwner(k+ndeltap)))*meshSf(1:3,k))
!
!				do i=1,3
!					delta=U(i,meshFaceOwner(k))-U(i,meshFaceOwner(k+ndeltap))
!					do j1=1,3
!						vtmp = DynVisc * meshSf(j1,k)*delta
!						do l=1,4
!							gradU(i,j1,meshFaces(l,k))=gradU(i,j1,meshFaces(l,k))+vtmp ! m2 * m/s
!							gradU(i,j1,meshFaces(l,k+ndeltap))=gradU(i,j1,meshFaces(l,k+ndeltap))-vtmp
!
!							if(i.eq.j)then
!								gradU(i,j1,meshFaces(l,k))=gradU(i,j1,meshFaces(l,k)) + BulkViscousStress ! m2 * m/s
!								gradU(i,j1,meshFaces(l,k+ndeltap))=gradU(i,j1,meshFaces(l,k+ndeltap)) - BulkViscousStress
!							end if
!						end do
!					end do
!				end do
!			end do
!		end if
!	end do
!
!	include 'netPointsInterface.h'
!
!	do k=1,nPoints
!		do i=1,3
!			do j=1,3
!				grad(i,j)=gradU(i,j,k)/meshVPoints(k) !  ! 1 / s
!			end do
!		end do
!		do i=1,3
!			do j=1,3
!				gradU(i,j,k) = grad(i,j)+grad(j,i) ! NOTE: include only the first viscosity, VSemiletov
!				!gradU(i,j,k)=nu*(grad(i,j)+grad(j,i)) ! NOTE: include only the first viscosity, VSemiletov
!					! [gradU(i,j,k)] = m2/s / s = m2/s2
!			end do
!		end do
!	end do
! ****************************************************************************************************************************** !
! viscosity 
!Counter = 0
! cycle on internal faces
!iCell = 0

!	sumViscousCellFluxes = 0d0

	do k=1,nInternalFaces
		vtmp=meshSf(1,k)*Usnew(1,k)+meshSf(2,k)*Usnew(2,k)+meshSf(3,k)*Usnew(3,k)
		i=meshFaceOwner(k)
		j=meshFaceNeighbour(k)
		!	write(*,*) 'internal k vtmp=',vtmp
		rhonew(i)=rhonew(i)+vtmp*Rhosnew(k)
		stress=0.0d0
		do l=1,3
			do m=1,3
				do m1=1,4
					stress(l,m) = stress(l,m) + gradU(l,m,meshFaces(m1,k)) / 4.0d0 ! kg/m3 * m2/s2 = kg / (m * s2)
!					stress(l,m)=stress(l,m)+Rhosnew(k)*gradU(l,m,meshFaces(m1,k))/4.0d0 ! kg/m3 * m2/s2 = kg / (m * s2)
				end do
			end do
		end do

!		faceGradU_view(1:3,1:3,k) = stress(1:3,1:3)

		DivVStress = 0.0d0
		do l = 1,3
			do m = 1,3
				DivVStress = DivVStress + meshSf(l,k) * Usnew(m,k) * stress(l,m)
			end do
		end do

		Enew(i)=Enew(i)+vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) + DivVStress
		Unew(1:3,i)=Unew(1:3,i)+vtmp*Rhosnew(k)*Usnew(1:3,k)+(Psnew(k)-1d5)*meshSf(1:3,k) &
			+meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)
		rhonew(j)=rhonew(j)-vtmp*Rhosnew(k)

		Enew(j)=Enew(j)-vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) - DivVStress
		Unew(1:3,j)=Unew(1:3,j)-vtmp*Rhosnew(k)*Usnew(1:3,k)-(Psnew(k)-1d5)*meshSf(1:3,k) & ! Psnew(k)*meshSf(1:3,k) &
			-meshSf(1,k)*stress(1,1:3)-meshSf(2,k)*stress(2,1:3)-meshSf(3,k)*stress(3,1:3)
				! [Unew] = m2 * stress = m2 * kg / (m * s2) = kg * m / s2

!		sumViscousCellFluxes(1:3,i) = sumViscousCellFluxes(1:3,i) - ( meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3))
!		sumViscousCellFluxes(1:3,j) = sumViscousCellFluxes(1:3,j) - ( -meshSf(1,k)*stress(1,1:3)-meshSf(2,k)*stress(2,1:3)-meshSf(3,k)*stress(3,1:3))

!		ViscousFluxes_view(1:3,1,k) = meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)
!		ViscousFluxes_view(1:3,2,k) = - meshSf(1,k)*stress(1,1:3)-meshSf(2,k)*stress(2,1:3)-meshSf(3,k)*stress(3,1:3)

!		if(nstep.eq.-1)then
!			if(abs(meshCellCenters(1,i)+0.05).lt.0.001.and. &
!				abs(meshCellCenters(2,i)+0.677).lt.0.001.and. &
!				abs(meshCellCenters(3,i)+0.05).lt.0.001)then
!
!				write(2,'(i3,i4,i4,8e13.5,3f8.3)') nstep,i,k, vtmp, Rhosnew(k), Usnew(1:3,k), Psnew(k), meshSf(2,k), &
!						meshFaceCenters(1:3,k)
!			end if
!			if(abs(meshCellCenters(1,j)+0.05).lt.0.001.and. &
!				abs(meshCellCenters(2,j)+0.677).lt.0.001.and. &
!				abs(meshCellCenters(3,j)+0.05).lt.0.001)then
!
!				write(2,'(i3,i4,i4,8e13.5,3f8.3)') nstep,j,k, -vtmp, Rhosnew(k), Usnew(1:3,k), Psnew(k), -meshSf(2,k), &
!						meshFaceCenters(1:3,k)
!			end if
!		end if

!		if(abs(meshCellCenters(1,i)+0.05).lt.0.001.and. &
!				abs(meshCellCenters(2,i)+0.677).lt.0.001.and. &
!				abs(meshCellCenters(3,i)+0.05).lt.0.001)then
!			write(*,*) '@@@@@ ',iCell, ' P*Fy=', Psnew(k)*meshSf(2,k)
!			iCell = iCell +1
!		end if	

!		if(abs(meshCellCenters(1,j)+0.05).lt.0.001.and. &
!					abs(meshCellCenters(2,j)+0.677).lt.0.001.and. &
!					abs(meshCellCenters(3,j)+0.05).lt.0.001)then
!				write(*,*) '@@@@@ ',iCell, ' P*Fy=', -Psnew(k)*meshSf(2,k)
!			iCell = iCell +1
!		end if
	end do		

! cycle on boundary faces
	do k=nInternalFaces+1,nFaces
		vtmp=meshSf(1,k)*Usnew(1,k)+meshSf(2,k)*Usnew(2,k)+meshSf(3,k)*Usnew(3,k)
		i=meshFaceOwner(k)
		stress=0.0d0
		do l=1,3
			do m=1,3
				do m1=1,4
					stress(l,m)=stress(l,m) + gradU(l,m,meshFaces(m1,k))/4.0d0
!					stress(l,m)=stress(l,m)+Rhosnew(k)*gradU(l,m,meshFaces(m1,k))/4.0d0
				end do
			end do
		end do

!		faceGradU_view(1:3,1:3,k) = stress(1:3,1:3)

		DivVStress = 0.0d0
		do l = 1,3
			do m = 1,3
				DivVStress = DivVStress + meshSf(l,k) * Usnew(m,k) * stress(l,m)
			end do
		end do

		rhonew(i)=rhonew(i)+vtmp*Rhosnew(k)
		Enew(i)=Enew(i)+vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) + DivVStress
		Unew(1:3,i)=Unew(1:3,i)+vtmp*Rhosnew(k)*Usnew(1:3,k)+ (Psnew(k)-1d5)*meshSf(1:3,k) & ! Psnew(k)*meshSf(1:3,k) &
				+meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)
		
!		sumViscousCellFluxes(1:3,i) = -(sumViscousCellFluxes(1:3,i) + meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3))
!		ViscousFluxes_view(1:3,1,k) = meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)

!		if(nstep.eq.-1)then
!			if(abs(meshCellCenters(1,i)+0.05).lt.0.001.and. &
!				abs(meshCellCenters(2,i)+0.677).lt.0.001.and. &
!				abs(meshCellCenters(3,i)+0.05).lt.0.001)then

!				write(2,'(i3,i4,i4,8e13.5,3f8.3)') nstep,i,k, &
!				vtmp, Rhosnew(k), Usnew(1:3,k), Psnew(k), meshSf(2,k), meshFaceCenters(1:3,k)
!			end if
!		end if

!		if(abs(meshCellCenters(1,i)+0.05).lt.0.001.and. &
!					abs(meshCellCenters(2,i)+0.677).lt.0.001.and. &
!					abs(meshCellCenters(3,i)+0.05).lt.0.001)then
!			write(*,*) '@@@@@ ',iCell, ' P*Fy=', Psnew(k)*meshSf(2,k)
!			iCell = iCell +1
!		end if
	end do	

!	do i = 1,nCells
!		if(abs(meshCellCenters(1,i)+0.05).lt.0.001.and. &
!				abs(meshCellCenters(2,i)+0.677).lt.0.001.and. &
!				abs(meshCellCenters(3,i)+0.05).lt.0.001)then
!				write(*,*) '@@@@@ Unew(1:3)=', Unew(1:3,i), ' Cell=', meshCellCenters(1:3,i)
!		end if
!	end do

	do i=1,nCells
		rhonew(i)=rho(i)-rhonew(i)/meshV(i)*dt2

! For testing acoustic post-processor !
!		memR = (meshCellCenters(1,i)+0.5d0)**2 + &
!			(meshCellCenters(2,i)+0.5d0)**2 + &
!			(meshCellCenters(3,i)+0.5d0)**2
!
!!		if(abs(meshCellCenters(1,i)+0.55d0).lt.0.001d0.and. &
!!				abs(meshCellCenters(2,i)+0.55d0).lt.0.001d0.and. &
!!				abs(meshCellCenters(3,i)+0.55d0).lt.0.001d0)then
!		if(memR.lt.ZeroEpsilon)then
!			rhonew(i) = rhonew(i) + & ! dt2 / meshV(i) * 0.001d0 * sin(50.0d0 * 3.1415926535897932384626433832795d0 * ttime)
!								  - 2d0 * 0.001d0/(50.0d0 * 3.1415926535897932384626433832795d0) * &
!								  sin(50.0d0 * 3.1415926535897932384626433832795d0 * (ttime+0.5d0*dt2)) * &
!								  sin(50.0d0 * 3.1415926535897932384626433832795d0 * 0.5d0*dt2)/ meshV(i)
!					It was really important integrating of RHS!
!		end if

		Enew(i) = (E(i)*rho(i)-Enew(i)/meshV(i)*dt2)/rhonew(i) ! E(i) !

		if(SourcePresent)then
			if(i.eq.iCellSource)then
				call CalcSource(ttime,Q)
				Q = Q*dt2  / meshV(i)
				Enew(i) = Enew(i) + Q / rhonew(i)
!				Rhonew(i) = Rhonew(i) + Q
!write(111,'(i4,4f25.10,24f15.5)') nstep,ttime, Q,Enew(i),P(i), &
!Rhosnew(meshCells(1:6,i)), Usnew(1,meshCells(1:6,i)), Psnew(meshCells(1:6,i)), Esnew(meshCells(1:6,i))
			end if
		end if

!		if(i.eq.1760)then ! 1875)then
!			Pi = 3.1415926535897932384626433832795
!		    ek = 1d7
!			a = 10.0d0 
!			Q = ek * sin(2.0d0 * a * Pi * ttime)
!			Enew(i) = Enew(i) + Q*dt2/rhonew(i)
!		end if
		Unew(1:3,i)=(U(1:3,i)*rho(i)-Unew(1:3,i)/meshV(i)*dt2)/rhonew(i)
			! [Unew] = kg * m / s2 / m3 * s / (kg/m3) = 1 / s2 / m2 * s * m3 = m / s

		pnew(i)=(gam-1.0d0)*(Enew(i)-(Unew(1,i)**2+Unew(2,i)**2+Unew(3,i)**2)/2.0d0)*rhonew(i)
		Tnew(i)=pnew(i)/(rhonew(i)*cv*(Gam-1.0d0))

		! determaine Currant's Time Step
		vtmp=DbCells(i)/(dsqrt(Unew(1,i)**2+Unew(2,i)**2+Unew(3,i)**2)+dsqrt(gam*pnew(i)/rhonew(i)))
		if (vtmp.lt.dtmin) dtmin=vtmp
	end do

!	if(SourcePresent)then
!		i = iCellSource
!		call CalcSource(ttime,Q)
!		Q = Q*dt2  / meshV(i)
!		Rhonew(i) = Rhonew(i) + Q
!		Rho0 = Rhonew(i)

!		pnew(i) = PFON + Q
!		Rhonew(i) = ROFON + Q/ (1.4*PFON/ROFON)
!		Unew(1,i) = 0d0
!		Unew(2,i) = 0d0
!		Unew(3,i) = 0d0
!		Enew(i) = pnew(i) / (0.4d0*Rhonew(i))
!		Tnew(i)=pnew(i)/(rhonew(i)*cv*(Gam-1.0d0))

!_______________________________________________________
!		Unew(1,i) = Unew(1,i) * Rho0 / (Rho0 - Q)
!		Unew(2,i) = Unew(2,i) * Rho0 / (Rho0 - Q)
!		Unew(3,i) = Unew(3,i) * Rho0 / (Rho0 - Q)

!		Enew(i) = Enew(i) * Rho0 / (Rho0 - Q)
!		Enew(i) = Enew(i) + Q
!_______________________________________________________
!		Unew(1,i) = Unew(1,i) + Q * U(1,i)
!		Unew(2,i) = Unew(2,i) + Q * U(2,i)
!		Unew(3,i) = Unew(3,i) + Q * U(3,i)
!		Enew(i) = Enew(i) + Q * E(i)
!	end if

!	if (ttime.gt.ttimeout) then
!		nout=nout+1
!	!	do j=1,8
!	!	xg(1:3)=xg(1:3)-meshPoints(1:3,meshCellPoints(j,l))/8.0d0
!	!	end do
!!		if (mod(nout,25).eq.0) then
!!			write(32) ttime,0,nCells,(Unew(1:3,i),pnew(i),Tnew(i),i=1,nCells)
!!			if (rank.eq.0) write(*,'(a18,3e12.5)') 'ttime,dtmin,dt2=',ttime,dtmin,dt2
!!		end if
!
!		do i=1,nCellsOut
!		!	vtmp=(rhonew(meshCellOut(i))-rho(meshCellOut(i)) )/dt2
!			vtmp=rhonew(meshCellOut(i))
!			Cellout(1:3,i)=Cellout(1:3,i)+Unew(1:3,meshCellOut(i))*(ttime-ttimeout)
!			Cellout(4,i)=Cellout(4,i)+pnew(meshCellOut(i))*(ttime-ttimeout)
!			Cellout(5,i)=Cellout(5,i)+vtmp*(ttime-ttimeout)
!		end do
!!	write(*,*) rank,nCellsOut
!		write(outFW) ttimeout,0,nCellsOut,((cellOut(j,i)/dtout,j=1,5),i=1,nCellsOut)
!		if (rank.eq.0) write(*,'(a18,i10,3e12.5)') 'nout,ttimeout=',nout,ttimeout
!		do i=1,nCellsOut
!		!	vtmp=(rhonew(meshCellOut(i))-rho(meshCellOut(i)) )/dt2
!			vtmp=rhonew(meshCellOut(i))
!			Cellout(1:3,i)=Unew(1:3,meshCellOut(i))*(dt2-(ttime-ttimeout))
!			Cellout(4,i)=pnew(meshCellOut(i))*(dt2-(ttime-ttimeout))
!			Cellout(5,i)=vtmp*(dt2-(ttime-ttimeout))
!		end do
!		ttimeout=ttimeout+dtOut
!	else
!		do i=1,nCellsOut
!		!	vtmp=(rhonew(meshCellOut(i))-rho(meshCellOut(i)) )/dt2
!			vtmp=rhonew(meshCellOut(i))
!			Cellout(1:3,i)=Cellout(1:3,i)+Unew(1:3,meshCellOut(i))*dt2
!			Cellout(4,i)=Cellout(4,i)+pnew(meshCellOut(i))*dt2
!			Cellout(5,i)=Cellout(5,i)+vtmp*dt2
!		end do
!	end if
